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We develop continuum field model for crack propagation in brittle amorphous solids. The model 
is represented by equations for elastic displacements combined with the order parameter equation 
which accounts for the dynamics of defects. This model captures all important phenomenology of 
crack propagation: crack initiation, propagation, dynamic fracture instability, sound emission, crack 
branching and fragmentation. 
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The dynamics of cracks is the long-standing chal- 
lenge in solid state physics and materials science 
The phenomenology of the crack propagation is well- 
established by recent experimental studies l^j-^: once 
a flux of energy to the crack tip passes the critical value, 
the crack becomes unstable, it begins to branch and emits 
sound. Although this rich phenomenology is consistent 
with the continuum theory, it fails to describe it because 
the way the macroscopic object breaks depends crucially 
on the details of cohesion on the microscopic scale . 

Significant progress in understanding of fracture dy- 
namics was made by large-scale (about 10^ atoms) molec- 
ular dynamics (MD) simulations |l^Jl^. Although lim- 
ited to sub-micron samples, these simulations were able 
to reproduce several key features of the crack propaga- 
tion, in particular, the initial acceleration of cracks and 
the onset of dynamic instability. However, detailed un- 
derstanding of the complex physics of the crack propaga- 
tion still remains a challenge 

The uniform motion of the crack is relatively well- 
understood in the framework of the continuum theory 
iQ. Most of the studies treat cracks as a front or in- 
terface separating broken/unbroken materials and propa- 
gating under the forces arising from elastic stresses in the 
bulk of material and additional cohesive stresses near the 
crack tip jTs]-^ . Although these investigations revealed 
some features of the oscillatory crack tip instability, they 
are based on built-in assumptions, e.g., on specific de- 
pendence of the fracture toughness on velocity, structure 
of the cohesive stress etc. To date there is no continuum 
model capable to describe in the same unified framework 
the whole phenomenology of the fractures, ranging from 
crack initiation to oscillations and branching. 

In this Letter we present a continuum field theory of 
the crack propagation. Our model is the wave equations 
for the elastic deformations combined with the equation 
for the order parameter, which is related to the concen- 
tration of material defects. The model captures all im- 
portant phenomenology: crack initiation by small pertur- 
bation, quasi-stationary propagation, instability of fast 
cracks, sound emission, branching and fragmentation. 

Model. Our model is a set of the elasto-dynamic equa- 



tions coupled to the equation for order parameter p. 
We define the order parameter as the relative concentra- 
tion of point defects in amorphous material (e.g., micro- 
voids). Outside the crack (no defects) p = 1 and p — 
inside the crack (all the atomic bonds are broken). At 
the crack surface p changes from to 1 on the scale much 
larger than the inter-atomic distance, justifying the con- 
tinuum description of the crack . 

We consider the two-dimensional geometry focusing on 
the so-called type-I crack mode, see Fig. ||. Equations of 
motion for an elastic media |20l are: 



da 



J = l,2. 
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Ui are the components of displacements, TyAu; accounts 
for viscous damping, -q is the viscosity coefficient [Q, po 
is the density of material. In the following we set p^ = 1. 



The stress tensor Oij is related to deformations via 
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where is the elastic strain tensor, E is the Young's 
modulus and a is the Poisson's ratio. To take into ac- 
count effect of plastic deformations inside the material 
we introduce dependence of E upon p. For simplicity we 
consider linear dependence E — Eop, where Eq is the 
regular Young's modulus. The term ~ I'p in Eq. (|^) 
accounts for the hydrostatic pressure created due to gen- 
eration of new defects, v is the constant which can be 
obtained from comparison with the experimental data. 
Although this term can be formally interpreted as the 
hydrostatic pressure due to thermal expansion, we stress 
that the thermal expansion may be not the only mecha- 
nism contributing to the magnitude of this term. There is 
experimental evidence that the temperature at the crack 
tip can be of the order of several hundred degrees [ p2[ . 
However, this would imply the concept of thermal equi- 
librium which is unlikely to be achieved at the crack tip. 

One can observe that Eqs. (|l|) are linear elasticity equa- 
tions for p = 1, i.e., outside the crack, and have trivial 
dynamics for p = (there is no dynamics inside crack) . 
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We assume that the order parameter p is governed 
by pure dissipative dynamics which can be derived from 
the "free-energy" type functional !F, i.e., p — —5T/5p. 
Following Landau ideas on phase transitions |2^, we 
adapt the simplest form for the "free energy" T ~ 
/ dxdy{D\\/ p\^ + 4>{p)), where the "local potential en- 
ergy" (j) has minima at p = and p = 1. Choosing the 
polynomial form for (j){p) we arrive at 

p^D^p^ ap{l - p)F{p, uii) + fip)^ui. (3) 

Coupling to the displacement field enters through the po- 
sition of the unstable fixed point defined by the function 
F{p, uii), where uu is the trace of the strain tensor. Con- 
straint imposed on F(p, uu) is such that it must have one 
zero in interval 1 > p > 0: F{pc,uu) = 0, 1 > pc > 0, 
and dpF{p = Pc,uu) < 0. The simplest form of F satisfy- 
ing this constraint is F{p, uu) = 1 — {b — puii)p, although 
any monotonic function of p will show the qualitatively 
similar behavior. Here /3 and p are material constants re- 
lated to such properties as crack toughness and strain to 
failure. Coefficients D and a can be set to 1 by rescaling 
t — > at, Xi — *■ vxi where = D /a. 

The last term in Eq. (^) represents the coupling of the 
order parameter to the velocity field ii and is responsible 
for the localized shrinkage of the crack due to material 
motion. Since the specific form of the function /(p) is 
irrelevant, we take /(p) = cp(f — p) to insure that / 
vanishes at p = (no material) and p = 1 (no defects), 
where c is the dimensionless material constant. From our 
simulations we have found that this term in Eq. ^ is 
crucial to maintain the sharp form of the crack tip. 
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FIG. 1. Schematic representation of fixed-grips loading 

Static solutions. Eqs. (|l|)-(^) have a dip-like solution 
corresponding to the open gap far behind the crack tip 
(a "groove" along x-axis for our geometry, see Fig.(l)). 
The static one-dimensional equations read: 



dy 



dy'^ 



dpUyy_^^ ^_^(1_P)(1_(5_^^^^)^)^0, (4) 



with the fixed-grips boundary conditions (BC): Uy{y = 
±L) = ±LS {5 — relative displacement), p{y = ±L) — 1, 
and dyp{y = 0) = 0. Exclusion of Uyy from Eqs. (Q) yields 



where C is a constant of integration (C dy/p{y) = 
LS), P ^ b/il + pC), and ^ = yy/l + pC. The solution 
to Eq.(0) satisfying the BC for L ^ oo is: 



y/iP + 1)(1 - /3/2) cosh(gVF^) + 2-f3 
v/(/3 + 1)(1 - (3/2) cosh(^^/^^) + 2(3-1' 



(6) 



This solution exists for I < /3 < 2. A deep and wide crack 
opening is attainable if 2— (3 = e <C f . In this case the BC 
for Uy can be reduced to SL = C{L + ^T^y3/{eb)) yielding 
an equation for e since C = {b/2 — l)/p. For e <C 1 the 
width of the crack opening d, defined as p(d/2) = 1/2, is 
d = ■y/2/61n(24/e). After exclusion of e one arrives at 



d = 




(7) 



The solution to Eq. (Q) exists only if 5 exceeds some 
critical value 6c given by 6c « {b/2~- l)/p, which leads to 
the relation between strain to failure 6c and the material 
parameters p and b. The logarithmic, instead of linear, 
dependence of crack opening on system size L in Eq. (|^) 
is a shortcoming of the model resulting from oversimpli- 
fied dependence of the function F on uu . 

To study the dynamics of cracks we perform numerical 
simulations with Eqs. (0)-(||). We use explicit second- 
order numerical scheme with the number of grid points 
up to 4000 X 800. Our model reproduces all important 
phenomenology: crack arrest below the critical stress, 
crack propagation above the critical stress, oscillations 
of the crack velocity, crack branching and fragmentation. 
Selected results are presented in Figs. 2-5. 

Quasi- stationary propagation. We considered crack 
propagation initiated from a long notch with the length 
of the order 100 units. At relatively small loadings 6 we 
have observed a quasi-stationary propagation (no oscilla- 
tions). The crack produces the stress concentration near 
the tip, and the stress is relaxed behind the tip, see Fig. 
||b. The distribution of shear (Figs. ||c,§) is close to that 
expected from the elasticity theory. The angular depen- 
dence Tixy{d) of the shear stress axy near the tip is close 
to the theoretical dependence (Jxy{r,9) ^ r^^^'^YP^y + 
where YP^y{9) = sin(6') cos(30/2) obtained for the infinite 
stationary crack. The discrepancy can be attributed to 
finite-size effects, velocity correction, etc. We computed 
the angle 9m of the maximum shear stress vs crack speed 
V normalized by the Rayleigh speed Vr. (see Fig. |c, 
^ . As one derives from the linear elasticity, the angle in- 
creases with the speed of the crack [Q , in an agreement 
with our numerical results. 



iyy^C/p, 9|p = p(l-p)(l-/3p), (5) 
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The calculated dependence of the crack tip velocity on 
the effective fracture energy G ~ (5^, shown on Fig. (||), 
demonstrates an excellent agreement with the experi- 
mental data from Ref. The instability of the crack 
occurs when the velocity becomes of the order of 55% 
of the Rayleigh speed for the parameters of Fig. ||a-e. 
For parameters of Fig. |^ we have found lower value of 
the critical velocity, namely about 32% of the Rayleigh 
speed. In all cases the instability manifests itself as pro- 
nounced velocity oscillations, crack branching and the 
sound emission from the crack tip. 

Our calculations indicate absence of the minimal crack 
velocity, the so-called velocity gap . The initial veloc- 
ity jump, seen experimentally as well as in some of our 
simulations (see Fig. ^) , is attributed to the fact that the 
initial crack (notch) is too short or too blunt. 
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FIG. 2. Color-coded images of p{x,y) (a); hydrostatic pres- 
sure p = -p{uxx + Uyy) (b); and shear piu^y + Uy^) (c). 
Domain size 800 x 200, number of grid points 1600 x 400, 
and 5 = 0.05, 77 = 0.25, Eo = 10, a = 0.2,6 = 2.25, c = 11, 
v = 2.3, ^ = 9.2; p(x,y) for unstable propagation at (5 = 0.089 
(d) and 5 = 0.11 (e). Domain size 1200 x 200, 2400 x 400 
grid points, (f) p{x,y) for propagation with fragmentation, 
Eq = 100,0- = 0.36, c = 16, = 2.3,^ = 54 and 5 = 0.03, 
Domain size 2000 x 400, 4000 x 800 grid points. 




FIG. 3. 6m vs crack velocity V for parameters of Fig. ^-c. 
Solid line shows 9m vs V from linear elasticity [Q, o show 
result of numerical simulations. Inset: angular dependence of 
shear stress T.xy (filled circles). The dependence from linear 
elasticity theory for infinite crack is shown in solid line. 
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FIG. 4. Crack velocity V/Vr vs G/Gc, where Gc s the onset 
value for crack propagation, o correspond to stable propaga- 
tion, X to unstable propagation, parameters the same as for 
Fig. l^-e, o show the experimental results from Ref. jsj nor- 
malized by Vr = 1300 m/s, the Rayleigh speed in PMMA in 
the frequency range of 0.6-1 MHz M. 

Instability of cracks propagation. Taking the suffi- 
ciently large values of i' and c and starting from short 
cracks with the large load we observed consecutive crack 
branching. Since Eqs. (|l|)-(||) are homogeneous, these 
secondary crack branches typically retract after the stress 
at the tip of the shorter crack relaxes. Although this re- 
tracting may indeed take place, e.g., in vacuum the small 
cracks may heal, the oxidation of the crack surface and 
lattice trapping would prevent cracks from healing. In 
order to model these effects, one can introduce an addi- 
tional field representing concentration of oxygen and then 
couple it to the order parameter. In some simulations, 
we multiplied r.h.s. of Eq. (^) by a monotonic function 
w{x — Xtip)'. w{x > 0) = 1 and w{x — )■ —00) 0, where 
xtip is the crack tip position. Thus, wc slowed down 
the evolution of p behind the crack tip, which, in turn, 
prevents secondary cracks from healing. We succeeded to 
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obtain realistic crack forms, see Fig. gd,e. For fast cracks 
the "freezing" is not a necessity, since the retraction is 
rather slow. Fig. ||f shows results without freezing: mas- 
sive crack branching along with crack healing are present. 

Far away from the crack tip we have registered oscil- 
lations of hydrostatic pressure (see Fig. ^, Inset), which 
is a clear indication of the sound emission by the crack 
tip. The sound waves reflected from boundaries may also 
induce velocity oscillations, but they do not provide a 
mechanism for branching p^ . Increase in the applied 
displacement S results in increase of amplitude and the 
number of secondary branches (cm. Figs. |d-f and|). 
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FIG. 5. The crack tip velocity V normalized by Rayleigh 
velocity Vr vs time for 5 = 0.089 (solid line) and S = 0.11 
(dashed line), parameters of Fig. ^,e. Inset: oscillations of 
pressure p far away from the crack tip for 5 — 0.089. 

Some estimates are in order. Our unit of length A is the 
width of the craze zone and is of the order of a micron 
in PMMA (l|]. The unit of time r is obtained from 
A ^ I^ and Rayleigh speed Vr ^ lO'^m/s, which gives 
T ^ \/Eo10~^ s. In experiments the characteristic 
time of velocity oscillations tq is of the order of 1 fis. Our 
model gives tq ~ lOV ~ 0.1 - 1 i^s for = 10 - 100. 

Conclusion. We have developed a continuum field the- 
ory of the crack propagation. The central element of our 
approach is the description of crack by the order parame- 
ter. The proposed approach enables us to avoid the stress 
singularity at the crack tip and to derive the tip instabil- 
ity. Our model is complimentary to MD simulations of 
cracks and allows for a description of fracture phenomena 
on large scales. The parameters of our model can be ob- 
tained from comparison with the experiment. It will be 
interesting to derive the order parameter equation from 
discrete models of crack propagation [^,^ . 
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